function [output] = SL_main(input_array)

% data input
method = input_array.method;
data = input_array.data;
alt_vec = input_array.alt_vec;

D = data.D;
X = data.X;
Y = data.Y;
switch lower(method)
    case 'ols'
        Z = D; % OLS is a special case of IV
    case 'iv'
        Z = data.Z;
end
n = length(D);
N = n/2;
R = ceil(N^(2/3));
b = ceil(N^(1/3));

% In IV cases, beta0_iv-beta0=(Z1'*M2*X1)^{-1}(Z1*M2*U)
M2 = eye(n)-X*((X'*X)\X');
S = (Z'*M2)'.*Y/(Z'*M2*D/n);
beta0_hat = mean(S); % IV estimate

ydep = S;
x_reg = ones(length(S),1);
        
p_null = SongLeungU(x_reg,ydep,R,b,1,0);

p_alt = 0*alt_vec;
for i_alt = 1 : length(alt_vec)
    ydep_alt = ydep+alt_vec(i_alt);
    p_alt(i_alt) = SongLeungU(x_reg,ydep_alt,R,b,1,0);
end

%% OUTPUT
output.estimate = beta0_hat;
output.p_null = p_null;
output.p_alt = p_alt;

